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Abstract 

We present a study of nonlinear localized excitations called discrete breathers in a 
superconducting array. These localized solutions were recently observed in Josephson- 
junction ladder arrays by two different experimental groups [1-3]. We review the 
experiments made by Trias et al. [1] We report the detection of different single-site 
and multi-site breather states and study the dynamics when changing the array bias 
current. By changing the temperature we can control the value of the damping (the 
Stewart-McCumber parameter) in the array, thus allowing an experimental study 
at different array parameters. We propose a simple DC circuit model to understand 
most of the features of the detected states. We have also compared this model and 
the experiments with simulations of the dynamics of the array. We show that the 
study of the resonances in the ladder and the use of harmonic balance techniques al- 
low for understanding of most of the numerical results. We have computed existence 
diagrams of breather solutions in our arrays, found resonant localized solutions and 
described the localized states in terms of vortex and antivortex motion. 
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1 Introduction 



Linear models of crystals have been instrumental in developing a physical 
understanding of the solid state. Thermodynamic properties such as specific 
heat, transport properties such as electron relaxation times or superconduc- 
tivity, and even interactions with radiation can be understood by modeling 
a crystal as a lattice of atoms with fixed harmonic coupling. This leads to 
the conventional phonon-like analysis with a basis of plane waves as normal 
modes. However, certain properties of solids, such as thermal expansion, can- 
not be understood in this linear model. For example, the elastic constants of 
the atomic interactions may depend on the temperature or the volume and so 
make the interaction non-linear. This is termed an anharmonic effect and the 
usual approach is to use a more generalized Taylor expansion for the lattice 
coupling that includes more than just the harmonic term. 

Until quite recently anharmonic effects were only studied as perturbations to 
the fully solvable harmonic model. Then it was discovered [4] and proved [5] 
that in classical Hamiltonian systems non-linearity may lead to localized vi- 
brations in the lattice that cannot be analyzed using the standard plane wave 
approach. These intrinsic localized modes are time periodic and spatially lo- 
calized solutions and have been termed discrete breathers (DB). Their am- 
plitudes oscillate around a few sites in the lattice and they do not depend 
on impurities for their localization. The study of DB has been extended to 
the dynamics of coupled rotor lattices [6] where the terms "rotating localized 
modes" or "rotobreathers" were introduced. DB have also been proven to ex- 
ist in the dynamics of dissipative systems [7]. In this case chaotic localized 
solutions have been discovered [8]. DB have been found to play an important 
role in conditions far from equilibrium [9] and recently have been studied in 
disordered systems [10,11]. Excellent reviews on the topic are [12] and [13]. 

Because DB are generic modes in many non-linear lattices, they are the ob- 
ject of great theoretical and numerical attention in many diverse fields like 
condensed matter physics [14-17], mechanical engineering [18-20] and bio- 
physics [21]. Only recently have the first experiments been performed which 
detect intrinsic localized modes in quasi-one-dimensional charge-density-wave 
compounds [22], antiferromagnetic anharmonic crystals [23] and supercon- 
ducting arrays [1-3]. 

The existence of DB in a Josephson-j unction (JJ) array was first proposed by 
Florfa et al. [15] in the study of the dynamics of an AC-biased anisotropic 

ladder array. Both, oscillating and rotating localized modes were simulated 
and studied in this system [24,25]. Later, rotobreathers were also numeri- 
cally studied in the dynamics of inductively coupled junctions in DC-biased 
arrays [26,27]. Arrays were then designed, fabricated and measured and DB 
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Fig. 1. Sketch of the IV curve of an underdamped Josephson junction, 
were found [1-3]. 

The Josephson junctions studied in this article are made from a superconducting 
insulator-superconductor (SIS) tunneling structure that because of the Joseph- 
son effect behave as solid-state non-hnear oscillators. In the framework of 
the resistively and capacitively shunted junction (RCSJ) model a single JJ is 
modeled by a parallel combination of an ideal junction, a capacitor C, and a 
resistance R. The current of the SIS tunnel junction is then 



(1) 



The ideal junction has a constitutive relation of Ij = IcSinip where (p is the 
gauge-invariant phase difference of the junction and the voltage V across a 
junction is F = {^o/2TT)dip/dt. Thus, the current of SIS junction, i, in reduced 
units is given by 

i = if + Tip + sin ip = M{p) . (2) 



This current is normalized by the junction's intrinsic critical current Jc and 
r, the damping, is usually referred to as the Stewart-McCumber parameter 
/3c = = 2t^IcCR^ I^Q ($0 is the flux quantum). Time is normahzed by 
T = ^J^qC /2t:Ic. The RCSJ model equation is isomorphic to the equation of 
a driven pendulum. The mass is normalized to one and the viscous damping 

is r. 

Fig. 1 sketches a typical current-voltage (IV) curve of a single junction. The 
junction is biased by a DC current and the average voltage is measured. The 

IV curve of a JJ presents two different physical states. The first is a zero- 
voltage state, the so-called superconducting or quiet state; a current flows 
through the junction but no voltage difference appears. This state exists for 
values of the current smaller than the junction critical current Ic- At this value 
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of the current the junction switches to the superconducting gap voltage, Vg. 
The gap voltage results from the breaking of Cooper-pairs and causes the 
junction's resistance, and thereby the damping, to change in a complicated 
nonlinear way. If we increase the applied current further, the junction reaches 
its normal state and it behaves as a resistor, These resistive states are 
also called rotating or whirling modes. As the current decreases the junction 
returns to the gap voltage and then to its zero- voltage state at the retrapping 
current, 7^. Thus, the resistive state occurs for values of the current above the 
junction retrapping current and coexists with the superconducting state for 
currents between the critical and the retrapping values. The amplitude of this 
hysteretic loop is governed by the value of the damping F. The behavior of 
the curve (Fig. 1) close to Vg can be modeled using the RCSJ model with an 
appropriate nonlinear voltage dependence for the resistance (see section 5.2) 
or by using other more sophisticated models [28]. 

We can design JJ arrays of different geometries and parameters. Networks of 
junctions are valuable model systems for the study of coupled non-linear oscil- 
lators. For instance, solid-state physical realizations of the Frenkel-Kontorova 
model for dislocations [29] and the two-dimensional X-Y model for phase tran- 
sitions in condensed matter [30] are two prominent examples. There has also 
been extensive studies on the soliton dynamics in one- dimensional Josephson 
arrays [31]. The experimental advantage of Josephson networks is that there is 
good control over the parameters because they are fabricated microelectronic 
solid-state circuits. Moreover, these networks can be designed for a wide range 
of oscillator parameters from the extremely underdamped to overdamped lim- 
its. 

Here we present an in-depth study of the experiments reported in [1]. We use 
more complete models to understand the dynamics of the array. In particu- 
lar we will not assume uniform current bias, the effects of temperature will 
be discussed, a nonlinear junction resistance is added, and a full-inductance 
matrix will be used. In the next section we introduce the governing equations 
and explain intrinsic localization in the ladder. In section 3 we will report on 
the experimental study of DB in our superconducting arrays. We develop in 
section 4 a simple circuit model which allows for the understanding of most of 
the experimental findings. Numerical simulations will be shown section 5. In 
section 6 we perform a linear analysis that yields resonance frequencies and 
decay lengths of excitations. Section 7.1 is devoted to a numerical study of 
single-site DB solutions in the ladder. There we compute regions of existence 
of DB at different array parameters, and we will analyze typical simulated 
IV curves for the cases of type-B and type-A DB. We also study the vortex 
patterns associated to both solutions. The major results are summarized in 
the final section. 
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Fig. 2. Anisotropic ladder array with uniform current injection. Vertical junctions 
(with superscript v) have critical current lev and horizontal junctions (with super- 
script t and b) have a critical current Ich- 



2 Josephson Ladders 



2.1 Ladder equations 



Fig. 2 shows the circuit diagram for the Josephson array in ladder geometry 
and with uniform current injection. The junctions are marked by an "x". 
Horizontal junctions have critical currents of Ich while vertical junctions have 
a critical current of lev Anisotropic arrays are fabricated by varying the area 
of the junctions. In our junctions, the critical current and capacitance are 
proportional to this area. Due to the constant IcRn product, the normal state 
resistance is inversely proportional to the junction area. The anisotropy pa- 
rameter h can then be defined as h — Ich/ lev — Ch/C^ — Ry/Rh- 

A ladder is a useful geometry to study DB because vertical junctions can 
play the role of pendula while the horizontal junction can act as controllable 
non-linear coupling. Thereby, a ladder can be roughly thought of as a one- 
dimensional chain of pendula that are coupled by non-linear springs. If the 
phases of the vertical junctions are interpreted as as particle coordinates, then 
the ladder is in essence a one-dimensional chain of particles with non-hnear on- 
site potential as well as non-linear interactions since the horizontal junctions 
arc nonlinear Josephson elements. It is known that lattices with non-linear on- 
site potential or non-linear interactions are needed to support DB excitations. 

The junctions in the array are coupled by means of current conservation and 
fluxoid quantization. Kirchhoff 's current conservation law (KCL) at the top 
node yields 

+ /ext - /] - /J = 0. (3) 



A consequence of the open boundary conditions is that the current on the top 
horizontal junctions must be equal but opposite to the current in the bottom 
horizontal junctions. Thus, /j" = — /* = Ij. We will normalize all the currents 



5 



by lev Also, we will refer to horizontal junctions by the superscript h (which 
is not to be confused with the anisotropy) when we are dealing with either 
the top or bottom horizontal junction. 

Fluxoid quantization causes the circulation of the gauge invariant phase dif- 
ferences around a loop to be equal to the flux of the total (external plus 
induced) magnetic field through the loop. When we impose this condition on 
one of our cells and assume only external and self-induced fields, that is, we 
neglect mutual inductances between different cells, we find 

- - ^5 + ^5 + 2^/ + - 0- (4) 

Here, / = $ext/^o is the normalized applied flux per unit cell and A = 
^q/2'kIcvLs where Lg is the self-inductance of the loop and if is the nor- 
malized mesh current so that if /\ is the normalized self-induced flux of a 
cell. 

The vorticity, rij is defined through the expression 

bl] - [V'l+i] - ¥,\ + [<^5] = 27r(n, - / - ff') (5) 

where [</?] represents the phases modulus 27r, and /j"'^ = if /2t:\. This expres- 
sion is equivalent to Eq. (4) and thus also referred to as fluxoid quantization. 

We let the functional N{^) = (p + + sin ^ represents the current through a 
junction in the RCSJ model. The resulting set of nonlinear coupled equations 
can be written as 

N{ip]) = \{ip)^, - 2^'. + ip)_, + if] - ip]_, - if] + + Zext 

= -IW] - - + ^) + 2^/} (6) 

We can identify in these equations a discrete Laplacian term V^v^J = V'J+i — 
2(/9j + which accounts for the interaction between vertical junctions, and 
discrete first order derivatives 5x^^j = 97^+1 — 9?} and (JxV^j-i = ^'j—^'j-i, which 
account for the interaction terms between vertical and horizontal junctions 
respectively. 

The external current is normalized as i^xt = hxt/hv The damping is F = 
^J(^Q/2nIcvR1Cv We note that because the anisotropy in our arrays is caused 
by varying the junction area, F is the same for every junction in the array. 
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In Eq. (6), j = 1 to and at the open boundaries, 



^Vn + 27r/ 

^l^^l-27:f (7) 

where the phases at j = and j = N + 1 are for mathematical convenience 
and do not represent real junctions. 

The physical currents through the junctions are /j = IchAf{^j) — —Ij — 
-Ichf^{ip'^ and I] = IcM{ip]). 

We can use these governing equations to compare the ladder with a JJ parallel 
array, a system with broad interest. In a parallel array the horizontal junc- 
tions of the ladder are replaced by superconducting wires which have a linear 
current-phase relation. In the ladder, the dynamics of the horizontal junc- 
tions can also be described by a linearized constitutive relation under some 
restricted circumstances. For instance, in the static case, when the horizon- 
tal junctions have no time dependence then — hsiiup^ h(p^. Prom the 
normalized Eq. (3), 

^J-W = ^5_l-^* (8) 



The right hand side is simply h{ip^j_i—ip^j) and we can find a similar relation for 
the bottom horizontal junctions. We can substitute these hnearized relations 
for (/?* and cp^ in Eq. (6) to get, 

•^(^I) = 7rT^^Vl + w (9) 



This is the discrete sinc-Gordon equation with a renormalized discreteness 
parameter of Aeg = h\/{h + 2\) and is equivalent to the governing equations 
of a Josephson- junction parallel array [31]. 

The difference between the parallel array and the ladder is obviously the exis- 
tence of horizontal junctions in the case of the ladder. The use of X^s to map 
the ladder to the discrete sine-Gordon model is correct only for the study of 
a reduced set of possible states of the array: Those for which only the convex 
part of the intcr-phase interaction is relevant. If we can neglect the dynamics 
of the horizontal junctions (for instance when studying static properties) then 
the above is a very good approximation as has been rigorously stated in [32]. 
The rotating DB states studied in this paper are a good example of a set of 
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Fig. 3. Different patterns for single-site DB solutions in the ladder. Each array is 
vertically biased by a constant DC current as shown in Fig. 2. 

solutions of the ladder which do not appear in the dynamics of the discrete 
sine- Gordon model. 



2.2 Localization in the Josephson ladder 

In Fig. 1 we showed that the IV curve of an underdamped JJ has an hysteresis 
loop for current values between the critical and the retrapping currents. In 
this current range the zero-voltage {V — 0) and rotating {V — Vg) attractors 
coexist, and it is this hysteresis loop that allows for the existence of breathers 
in the ladder with DC bias current. In the full ladder, the phase space is much 
more complex. However, when the vertical junctions are weakly coupled then 
it is possible that each of the junction attractors is essentially independent. 
A breather is then a localized state where one vertical junction is rotating 
while the others are in the zero voltage state. Under general considerations of 
nonlinear oscillators, it is not obvious that the phase space of the coupled 
system supports attractors of localized solutions. Indeed, in linear systems 
even when the oscillators are weakly coupled, the phase space does not support 
stable localized solutions. 

The coupling of vertical junctions occurs through horizontal junctions, geo- 
metrical inductances, and fluxoid quantization. But the most important con- 
tribution is from horizontal junctions and their influence is measured by the 
parameter h. If the couphng is too strong, i.e. h is too large, then localized 
solution cannot exist. We numerically found that small A and h — 0.25 are 
adequate values of the coupling parameters for studying DB [27]. 

Fig. 3 shows some of the possible conflgurations for single-site breathers in our 
ladders. The arrows indicate voltage polarity. Junctions without arrows are in 
the zero- voltage state. The solutions in Fig. 3 are single-site breathers because 
only one vertical junction is rotating. We also see that due to Kirchhoff's 
voltage law (KVL), there must be at least one other junction rotating in each 
of the neighboring cells. 
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Fig. 4. Simulation of 9x1 array with A = 0.05, F = 0.1 and h = 0.25. We have 
plotted the absolute value of the DC flux per unit cell at I = 0.7. The flux decays 
exponentially with a decay length of 0.32 for both solutions. 



The actual number and pattern of rotating horizontal junctions will determine 
the type of breather. We call breathers that have two rotating horizontal 
junctions, like Fig. 3(a) and Fig. 3(c), type A breathers. We see that in this 
case the voltages of all the rotating junctions are the same and so the breather 
solution is highly constrained. We will use the term asymmetric to refer to this 
set of type A solutions. 

In the case of the type B solution, the four nearest horizontal junctions rotate 
as depicted in Fig. 3(b). We will use the term symmetric to refer to type 
B solution. However, it is important to note that symmetric refers only to 
the fact that top and bottom horizontal junction in the same cell rotate. It 
does not always correspond to true up-down symmetric solutions. An up-down 
symmetric solution is a solution for which <^*(t) = — (/?^(t) mod27r. This is a 
solution of the system as can be inferred from Eqs. (6), where it is clear 
that A/"(93*) = —N'{ip'j). However, for certain values of the parameters, the 
underdamped character of the junctions allows for different solutions, such as 
the type A breathers, that do not obey the up-down symmetry condition. Also, 
for most type B solutions the magnitude of the voltage of a horizontal rotating 
junction is one half the voltage of the rotating vertical junction. Again, not 
all type B solutions obey this condition. 

Single-site breather solutions that have 3 horizontal junctions rotating such as 
those shown in Fig. 3(d) arc called hybrids of type A and B. The other possible 
single-site breathers that are not shown can also be classified as either A, B 
or hybrid. 

In the rest of this paper we will designate librating junctions as shorts and 
junctions that are rotating by an x as shown at the top of Fig. 4. Each branch 
represents a junction, but only the ones with a cross are rotating. Fig. 4(a) 
shows a type A breather and a plot of the numerically calculated DC flux per 
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Fig. 5. Schematic of the DB found experimentally in Fig. 8. Branches with x's 
depict rotating junctions. Bl through B8 represent type B breathers with the indi- 
cated number of rotating vertical junctions. State W is the whirling state were all 
of the vertical junction are rotating. 

unit cells for a 9 junction ladder array. Fig. 4(b) shows the decay of the DC 
flux for a type B breather. The flux decays exponentially so that our breather 
solutions are exponentially localized. In the case of the hybrid breather shown 
in Fig. 3(d), the decay of the flux in the right side of the array is the same 
as for the type B shown in Fig. 4(b) and the decay on the left side of the 
array is the same as the type A shown in Fig. 4(a). Note that we have plotted 
the modulus of the flux. The fluxes for all the cells of one side of the vertical 
rotating junction have opposite sign to the fluxes of the cells in the other side. 



Fig. 5 shows different type B breather solutions for which a set of contiguous 
vertical junctions are in the rotating state. We call them multi-site or m- 

site solutions where m refers to the number of vertical junctions rotating in 
the array. In general many different solutions are allowed where each vertical 
junction can be in the resistive or in the superconducting state while one or 
both horizontal junctions between vertical ones in different states also rotate. 
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Fig. 6. Current-voltage characteristic of an anisotropic Josephson junction ar- 
ray when no breathers are excited. The hysteresis between the depinning current 
(-^dep ^ 2 mA) and the retrapping current (/^ ~ 0.2 mA) is shown. Inset: Schematic 
of the anisotropic ladder array showing the bias circuit. Iq is the total applied cur- 
rent and Rb the bias resistances are 250. 

The above discussion focused on the existence of rotobreathers in DC-biased 
arrays. In the case of AC-biased arrays the JJ ladder in addition to rotating 
localized modes also supports oscillating localized modes or oscillobreathers 
where one vertical junction describes a large amplitude oscillation when com- 
paring with other vertical junctions. Such modes were studied in [15,24,25] 
for non-inductive arrays. These modes persist when inductances are added to 
the model and also different type-B and type-A families of AC-biased roto- 
breathers can be identified [33] . 



3 Experiments 

3.1 Experimental observation of breathers 

We have designed and measured several anisotropic ladders. The inset of Fig. 6 
shows a schematic of the measured arrays. The junctions are fabricated using 
a Nb-Al20a;-Nb tri-layer technology with a critical current density of about 
IkA/cm^ [34]. The current is injected and extracted through bias resistors in 
order to distribute it as uniformly as possible in the array. These resistors are 
large enough so as to minimize any deleterious effects on the dynamics. Our 
ladders have 3x3 /xm^ horizontal junctions and 6x6 yum^ vertical junctions. 
Vertical junctions have been designed with four times the area of the horizontal 
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ones. Thus, the anisotropy ratio h is approximately 0.25. The bias resistors 
are 25 Vt. 

There are voltage probes in the fourth, fifth, sixth and ninth vertical junctions 
to measure V4, V^, Vq and V^. The voltage probes can also be used to measure 
the top horizontal junctions in the middle, which we denote as V^t and V5T, 
or any other combination of terminals. 

Prom the measured normal state resistance we calculate lev = 360 /^A and 
Ich — 90 at T = K. The dimensionless penetration depth A, which mea- 
sures the inductive coupling in the array, is defined as ^o/2tiLsIcv We estimate 
the mesh inductance Lg — 30.2 pH from numerical modeling [35] so A = 0.04 
at T = K. 

To determine F we need to measure the subgap resistance. Different ap- 
proaches can be used. One possibility is to calculate this damping from the 
measured return current of the junction. The model used for the return current 
determines the subgap dissipation and resistance [36-39]. We will estimate 
the value of V in our experiments from the retrapping current by the re- 
lation Ir/NIcv — (4/7r)r, where is the number of vertical junctions. This 
expression can be found from a simple energy argument valid for small values 
of the damping. The energy injected into the system by the applied current 
must equal the energy to "rotate" the junction one full period [28]. Thus from 
the experiment shown in Fig. 6 we infer a value of F ~ 0.08 {(3c ~ 160) at 
T = 5.2K. 

We show a typical IV curve of a ladder without DB in Fig. 6. We measure 
the time- averaged voltage of the 9-th vertical junction as a function of the 
uniformly applied current. The junction is in a zero- voltage state as we increase 
the applied current from zero. When the applied current reaches the depinning 
current /^ep at about 2 mA the junction switches from zero- voltage state to 
the superconducting gap voltage, Vg, which at this temperature is 2.5 mV. If 
we increase the applied current further, the junction reaches its normal state 
and it behaves as a resistor, of 5^2. As the current decreases the junction 
returns to the gap voltage and then to its zero- voltage state at the retrapping 
current, ~ 0.2 mA. 

Sometimes, when we sweep the applied current we find that DB solutions 
appear spontaneously: They can be thermally excited when the applied current 
is close to /dep- However for our experiments, we have developed a simple 
reproducible method of exciting a breather: (i) Bias the array uniformly to 
a current below depinning current; (ii) increase the current injected into the 
middle vertical junction (V5) until its voltage switches to the gap; (iii) reduce 
this extra current in the middle junction to zero. Other procedures are possible. 
For instance, we can increase the current for the middle vertical junction first 
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Fig. 7. Measurement of the time-avcragcd voltages of five junctions of the array with 
the breather as the applied current is increased at T = 5.2 K. We first biased the 
ladder at 1.4 mA and excited a breather as indicated in the text. Then, the applied 
current is increased. Below /+ ^ 1.95 mA we see the signature measurement of the 
breather and above the breather becomes unstable and the array switches to the 
whirling state. 

until it rotates and then increase the array bias current as was described in [2]. 
We could also inject this extra current into a horizontal junction. All of these 
methods produce breathers in our arrays and thereby hinting at the generic 
nature of breather solutions in our ladders 

Fig. 7 shows the result after we have excited the breather and we have in- 
creased the array current. The breather is excited at ~ 1.4 mA and then 
the junction voltages are measured as the applied current is increased. We 
simultaneously measure the voltages of the vertical junctions (V4, V5 and Vq) 
and the top two horizontal junctions, V^t and V^t- The DB is localized in 
the fifth vertical junction and is a type B breather since the top, V^t, and 
bottom, V5B, horizontal junctions have opposite voltage. We also find that the 
horizontal voltages are half in magnitude to the rotating vertical junction V5. 
This is the Bl breather state shown in Fig. 5. 

The breather exists until a maximum current /_|_ ~ 1.95 mA is reached. If the 
applied current is further increased then Fig. 7 shows that the voltages of 
the 4th and 6th vertical junctions jump to the gap voltage while those of the 
horizontal junctions go to zero. In fact all the vertical junction have jumped to 
the gap voltage. The ladder has now all of the vertical junctions rotating and it 
is in the "whirling" state. The state is depicted as W in Fig. 5. In section 4 we 
will present a circuit model that will relate /+ to the array depinning current. 

In Fig. 8 we measure the voltage as the apphed external current is decreased. 
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Fig. 8. Measurement of time- averaged voltages as the applied current is decreased 
at T = 6.0 K. We show the voltages of three vertical junctions (V4,V5 and Vq) and 
the voltage measured in the top branch between the middle and one of the edges 
of the array (Vtop)- We first biased the ladder at 1.25 mA and excited a breather 
as indicated in the text. Then, the applied current is decreased. The nine steps 
corresponds to different type B m-site breathers. The dashed lines are the expected 
minimum currents based on a retrapping model [Eq. (17)]. The inset shows the 
values used to fit the minimum currents. 



We excite the breather at la ~ 1.25 mA and then the voltages are measured 
as the current is decreased. We have measured vertical junctions four, five 
and six. We have also measured the total voltage of the horizontal junctions 
five through eight and this sum voltage is referred as V^op- This allows us to 
reconstruct the m-site breather state in terms of the rotating junctions of the 
array. 

Our experiments show that as we decrease the applied current the single-site 
breather will usually decay into an m-site breather state. From /„ 1.25 mA 
to la ~ 1.05 mA the vertical voltages V4 = Vq = 0. This is the single-site type B 
breather, the state depicted as Bl in Fig. 5. If we further decrease the applied 
current from 7_ = 1.05 mA, V4 jumps from zero to the gap voltage. This is a 
two site breather shown schematically in Fig. 5 as B2. As we further decrease 
the current we can count nine discontinuous curves, each one corresponding 
to the switching of a vertical junction. At 0.3 mA all of the vertical junctions 
return to their zero-voltage state via a retrapping mechanism analogous to 
that of a single pendulum. 

From these experiments wc conclude that this shifting of the voltage corre- 
sponds to at least one vertical junction switching from the zero- voltage state 
to the rotating state. If we assume that only consecutive junctions switch, then 
every curve in the measurement can be associated with one m-site breather. 
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Fig. 9. IV showing type Bl and type A breather at 4.9 K. Dashed Une at la ~ 2.3 mA 
is the maximum current for the type B breather and ~ 0.7 mA is the maximum 
current for the type A breather. 

Since we have measured Vtop we can reconstruct the ladder solutions. Fig. 5 
shows a schematic of the states measured in Fig. 8. 

The shapes of the IV curves in this multi-site breather regime are influenced by 
the junction nonlinear resistance and the redistribution of current when each 
vertical junction switches. This redistribution may also govern the evolution of 
the system after each transition to one of the other possible breather attractors 
in the phase space of the array. 

Fig. 8 shows a measurement at T = 6. OK and we find 10 different states can be 
distinguished (8 m-site breather states, the whirling state, and the zero voltage 
state). In general we only see three or four m-site breathers as we decrease 
the current, and often we see different ones even under similar experimental 
conditions. 

Fig. 9 shows an IV which includes a type A breather. First, a type B breather 
is excited at la ~ 0.9 mA and the current is decreased. As we decrease the 
current V4 and Vq have zero voltage while V5 = 2V57'. At ~ 0.58 mA the 
type B breather becomes unstable. In most of our measurements, the array 
will then suddenly switch to an m site breather as shown in Fig. 8, but in 
this case we find a new switching behavior: The type B breather switches 
to a type A solution. This state is a type A breather because the voltage 
^5 = ^5T and V4 and Vg are zero. As the current decreases the type A breather 
disappears at la ~ 0.53mA and the array is in the superconducting state. 
In our experiments, the ladder always returns to the superconducting state 
whenever a type A breather reaches its minimum current. 
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Fig. 10. /+ and /_ for Bl breather as a function of the applied magnetic field at 
r = 6.0K. 

As the current decreases the type A breather is only accessible for a small 
current range of 0.05 mA. However, it is possible to bias on the type A breather 
and increase the current to trace out the hysteresis loop. The tracing of the 
type A breather voltage step is also shown in Fig. 9. We see that the type 
A breather exists up to a current of ^ 0.7 mA. Once it becomes unstable the 
array dynamics usually jumps back to the Bl breather. 



3.2 Temperature and magnetic field dependence 



By sweeping the temperature and magnetic field we can study how the cur- 
rent range in which our breather exists is affected by a change of the array 
parameters. We define and study the evolution of four current values of impor- 
tance: the current when the array returns to the zero- voltage, Ir] the maximum 
zero- voltage state current, Jdep; and the maximum and minimum current for 
a breather state, /+ and /_. 

Fig. 10 shows the typical dependence of /+ and /_ for a Bl breather as a 
function of applied magnetic field. The IV curves were measured by applying 
a perpendicular magnetic field of to 300 mG using a magnetic coil that is 
mounted on the radiation shield of our probe. There is some / dependence. 
This can be expected since /+ is related to the array depinning current. But 
at this temperature A = 0.06 so we expect to have a large Meissner current 
and a correspondingly relatively fiat Jc dependence vs /. If A were larger 
then the breather dynamics and the Ic should show a stronger magnetic field 
dependence. 
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Fig. 11. Maximum (solid squares) and minimum (solid circles) currents for type 
B breathers, and maximum (open squares) and minimum (open circles) currents 
for type A breathers. Triangles are the ladder retrapping current. Lines through 
the solid nd open circles are the expected minimum current /_ from Eq. (17). The 
dashed line is the expected uncorrected maximum current from Eq. (18) while the 
solid line above it is the corrected value. 



We can study further the existence region of the breather by changing the 
temperature of the sample. In this way, we can vary the parameters to a 
certain degree. The temperature causes the lev of the junction to change and 
hence change F and A. For our arrays, the junction parameters can range from 
0.031 < F < 0.61 and 0.04 < A < 0.43 as the temperature varies from 4.2 K 
to 9.2 K. 

Fig. 11 shows how the maximum and minimum current of both type A and 
type B breathers are affected as we vary F. In Fig. 11 F < 0.2 corresponds to 
T < 6.7 K and A ~ 0.05. At these low temperatures, lev essentially remains 
constant so A docs not vary. However, the subgap resistance varies substan- 
tially as can be seen from the retrapping current measurements in Fig. 11. 
Therefore, there is a larger variation in F. Figure also shows a nice agreement 
in between the experiments (points) and the theory (lines) except for the case 
of the maximum current of the type A solutions. Also, experimentally we did 
not find breathers for F > 0.2. 



4 Circuit model 
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4-1 Introduction 



In this section, we will develop a simplified DC circuit model of our array in 
order to understand the region of existence of the breather solutions in our 
experiments. Also, this model will allow for evaluating the effect of the bias 
resistors in the dynamics of the array. 

In the experiment we apply the external current through bias resistors as 
shown in the inset of Fig. 6 in order to distribute it uniformly. If the resistors 
are very small, the horizontal junctions are effectively shunted by a small 
resistance and no DB solutions can be excited. If the resistors are large so 
that they dominate over every other impedance, then applied current will be 
almost uniform throughout the array. However, one drawback of using large 
bias resistors is that they will create local heating of the sample and affect the 
measurements. 

Thus, when studying DB, there are at least two complications with the bias 
resistors. The first is that when we excite a breather state, the applied current 
is not completely uniform. So questions arise about how this non-uniformity 
affects localization. Also as we decrease the current we see transitions between 
different m-site breather states. In each of these transitions vertical junctions 
switch from the zero voltage state to the rotating state. This switching causes a 
change in the junction impedance and thereby affects the current distribution. 
This redistribution might be important to understand the dominant drive of 
the pattern selection process between m-site breathers. However, we will show 
below that the effect of the bias resistors only adds a small correction to the 
calculation of the existence region of breather states. 

To get some physical intuition we will use a simple model where rotating 
vertical junctions have a resistance of and rotating horizontal junctions 
have a resistance of Rh- Librating junction will be modeled as shorts. We will 
reduce the array to a simple network of resistors and calculate DC properties. 
The equivalent resistor network for a single-site symmetric breather located 
on junction 5 in our 9 junction array is shown in Fig. 12. 

When the array is in the superconducting or in the whirling states, the applied 
current distributes uniformly across the vertical junctions. This will not be the 
case when we have a breather since some junctions are in the resistive state 
while others are superconducting. Moreover, when we have a breather and we 
measure the voltage of the fifth vertical junction (Fig. 8), we find that its volt- 
age shifts to a higher value when new vertical junctions switch to the rotating 
state. These shifts are identified as jumps in the effective current biasing the 
junction due to the redistribution of currents in the array. Roughly speak- 
ing, most of the applied current wants to fiow through the superconducting 
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Fig. 12. Equivalent DC circuit for single-site type B breather located at junction 5 
in a 9 junction array. Nodes are labeled a, b, and c. 

junctions. Whenever a vertical junction switches from a superconducting to a 
resistive state there is some extra current that is distributed throughout the 
array and thereby the effective bias of the fifth junction becomes a little larger. 
This extra bias results in a jump of the measured voltage to a larger value. 
The size of the voltage jump is dependent on the non-linear subgap resistance 
and is usually large even for small changes of the effective bias current. 



4-2 IV curves and current distribution 



In this section we will focus on interpreting this redistribution and shifting of 
the voltage by using a simple DC model of the ladder. Wc use to represent 
the total applied current in the array. This is the current of the experimental 
current source. The current applied to a particular junction through the bias 
resistors, i.e. the current through the bias resistor, will be designated Ij. We 
will use Ii to denote the current from node a to node b. This current Ii is the 
sum of the currents Ji through J4. Also, because all of our bias resistors are 
the same and due to KVL Ii = I2 = I3 = I^. So, // = 4/4 for instance. The 
current through node a to node c is since the breather is located in junction 
5. 

For the type B breather we measured, the voltage is twice as large as the 
voltage of the horizontal junction 

2Vn = K- (10) 



Therefore 



h = 2|^7,, (11) 
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and we will assume that = hRh. 



The circuit has left-right symmetry so the current through the right branches 
must equal the current through the left. KCL at node a yields la = ^Ii + J5 
where la is the total current applied, and at node c yields = 21 ^ + Iv KVL 
on the top left mesh gives IiRb/i — IhRh — hRh = 0. Combining the equations 
results in 



To generalize to type A breathers we just note that the horizontal junction 
voltage is the same as the vertical junction voltage Vh = V^. We can write 
sVh = Vj, where s = 1 for type A breathers and s = 2 for type B breathers. 
We can also generalize for m-site breathers by modifying the equivalent DC 
circuit accordingly. For instance, when m = 2 both rotating vertical junctions 
can be lumped into an equivalent impedance of Ry/2. Generalizing our circuit 
results in 



We can also calculate the IV curve of a vertical rotating junction by substi- 
tuting Ifi = {h/s)I^ = {h/s)V^/Rv in Eq. (13). The result is 

^=|l + ^+fl-I^)^|^. (14) 
N \ sm \ N ) sRb ] Rv 



Another important variable is the amount of effective current biasing a vertical 
rotating junction, say /s. 



4 

N 




(15) 



From our experiments Rh/Rb ~ 0.8, h = 0.25 and = 9. Thus for the case 
of a type-B single-site breather (m = 1, s = 2) we get J5 = 0.934/^/9 and 
Ij = 1.008/a/9 (j 7^ 5). Thus the vertical rotating junction is biased by a 
smaller DC current than the quiet vertical junctions. Also this non uniformity 
disappears as the bias resistance is made larger. 

^ These equations are valid for type A and type B solutions with m consecutive 
rotating vertical junctions and with junction 1 and N in the no rotating state. In 
a similar way, we can compute the equations for the case when junction 1 or 
is rotating, for hybrid breathers, or for multi-site breathers for which the rotating 
junctions are not consecutive. 
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4-3 Estimation of I- 



To calculate we assume that the instabihty that destroys an m-site breather 
is due to a junction retrapping mechanism. Even if it appears that some other 
mechanism is important, like a resonance, this instability occurs in the subgap 
region of the junction where the voltage varies rapidly while the current re- 
mains relatively constant. Therefore, a retrapping current should give a good 
estimate of /_ regardless of the physical mechanism. 

Prom simple energy consideration of an isolated junction, the retrapping cur- 
rent can be estimated as ATIc/tt. Then, when the horizontal junction reaches 
its retrapping current Ih = 4rich/'n'. Thus, ly = ATIchs/nh = iTIcvs/Ti, and 
the vertical junction is at s times the junction retrapping current. Conversely, 
if the vertical junction is at the retrapping current, then the horizontal junction 
is at 1/s of the junction retrapping current. Therefore, as the applied current 
decreases, the horizontal junction always reaches the retrapping current first 
for a type B breather and both the horizontal and vertical junctions reach the 
retrapping at the same time for a type A breather. This is the reason why type 
B breathers can decay into m-site breathers while type A breathers apparently 
decay into the superconducting state. When a type A breather reaches /_ all 
of the rotating junctions rctrap and the resulting state is more likely to be the 
superconducting state. On the other hand, when a type B breather reaches 
/_ only the horizontal junction retraps while the vertical junctions remain 
whirling. 

Prom Eq. (13), 




when Ih reaches its retrapping value. 

To use this formula we need the horizontal junction parameters. We can esti- 
mate Ch from the specific capacitance of the tri-layer and the junction area. We 
find Ch — 300 fP. Prom the constant IcRn product, we find that Ich = 90 /xA. 
However, the remaining parameter Rh is more difficult to estimate since it de- 
pends highly on the nonlinear subgap region. Instead of trying to calculate Rh 
directly, we will fit Eq. (16) to our measured /_ using Rh as our fitting parame- 
ter. We include the effect of Rh through the definition of F = ^J'^^JOmlch^fih- 
Then from Eq. (16), 

^^^^ 



21 



and Rh as the only free parameter. This is how we estimate Rh, or equivalently 
r, for a given temperature. 

The inset of Fig. 8 shows the fitted values for R^ and the dashed hues in Fig. 8 
show the resulting 7_. We note that for solutions B7, B8, and W we used a 
different equivalent circuit based on the schematic of Fig. 5. Wc sec that Rh is 
approximately 50 Q for all m. This value is not totally unexpected. Roughly 
speaking, the horizontal junctions are shunted by two bias resistor of 25 fl 
each and the retrapping current depends strongly on the equivalent junction 
impedance. Since the subgap resistance for our junctions can be several kft 
the equivalent horizontal junction impedance will be dominated by the shunts 
which add up to 50 fl. 

We can also easily understand why 7_ decreases as m increases. For a constant 
array bias at m = 1 some fraction of the current will flow through both 
horizontal junctions. At the same bias current for m = 2, there is a horizontal 
junction that is not rotating in between the rotating horizontal junctions. 
By symmetry considerations there is no current flowing through this quiet 
horizontal junction. Since the applied bias current is the same, this implies 
that a larger fraction of the applied current must flow through the rotating 
horizontal junctions when m = 2 than when m — 1. Thereby, /_ is smaller for 
m = 2. 



4-4 Estimation of 1+ 

To calculate 7+ we look for the librating junction that supports the maximum 
current. When this junction reaches Ic, the breather has reached its maximum 
current. It is straightforward to find that the critical junction is the first ver- 
tical junction that is not rotating (the one nearest to the rotating ones). Let 
/* be the current through the junction. KCL at node b of the array yields 
7* = Ih + h- Here we have assumed that there is not current in the horizontal 
quiet junctions. Since Ii = AI^ we can substitute for the currents to solve for 
/* in terms of /„. We can also generalize to m-site breathers. The result is 



The maximum applied current 7+ occurs when /* = lev In the limit Rb ^ Rh, 



la _ 2h/m + s + h{l- m/N) Rh/Rb j^ 
N ~ h{l + 2/m) + S + hRh/Rb 



(18) 



/+ 2h + sm 



I. 



(19) 



N h{m + 2) + sm 



cv 
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This current will underestimate the actual value for /_|_. This is because we 
have not taken into account any of the horizontal junction currents that are 
in the quiet state. 

The effect of fluxoid quantization is to redistribute the currents of the quiet 

junctions. The currents in the bias resistors and the rotating junctions will 
remain unaffected. If we consider the effect of the next nearest mesh to the 
breather, then KCL at node b will be Ia+Ih = I*+Ich sin(/c/i). Here, Ich sin(A;/j) 
is the current in the next horizontal junction. Fluxoid quantization in this quiet 
mesh yields k* — ky — 2kh — when / = and we neglect the self induced field. 
Here k* = sin^^ {I* / lev) and k^ is the phase of the next quite vertical junction. 
Adding self fields will tend to decrease Ich sin(A;/i) because smaller screening 
currents are needed as the inductance becomes larger. This correction will 
then tend to overestimate /+. We note that to calculate /+ we set /* = lev so 
k* = 7r/2. 

With fluxoid quantization kh = 7r/4 — A;^/2 at la = /+. To first order we expect 
the current in the quiet vertical junction to be la/N and kv = sin^^ {la/ NI^v)- 
For consistency, we apply KCL so I3 + IchSin{kh) — IcvSin{kv)- Again, we 
neglect the current of the next quiet horizontal junctions. We can solve for 
ky in terms of kh and substitute back into the fluxoid quantization condition. 
Using sinx X and cos a; 1 — x^/2 we can get a closed expression 

_ -h + ^h' + 8{l-h/Icv) 



This expression is only valid when /„ « /. 



+ • 



The current I3 equals I4 — Ii/4 and can be calculated from Fig. 12 and is 
simply 

2/m + s/h + Rh/Rb la ^21) 



2/m + s/h + (1 - m/N)Rh/RbN' 



In the limit Rb Rh h ^ h/N. 

The maximum current will now be Eq. (18) when /* = I^v + /c/tSin(A;ft) with 
kh defined in Eq. (20) and Eq. (21). The resulting equation is transcendental. 
But we know the correction due to kh should be small so we can linearize the 
sine term. The equation then has dependencies of on both sides but can be 
solved by isolating the square root and squaring. This results in a quadratic 
equation in with easily extractable roots. The equation can also be solved 
iteratively. One iteration usually results in a good approximation. Hence, if /° 
is the uncorrected value from Eq. (18) when 7* = lev, then the first correction 
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can be calculated when I* — lev + Ich sin(/c^) with = 7° substituted in 
in Eq. (21). 

4-5 Comparison with the experiments 

The result of the simplified DC circuit model is compared to the measurements 
in Fig. 8 and Fig. 11. We used the measurement in Fig. 8 to estimate the values 
of Rh and we use Eq. (17) to calculate /_ for the type A and B breathers. In 
Fig. 11, we see that the result agrees quite well for I . We have also calculated 
/+ from our uncorrected expression Eq. (18) and with the correction due 
to Meissner currents. The measured and expected results agree for type B 
breathers. However, the measured /+ for the type A breathers is much lower 
than what our circuit model predicts. 

We also note that the efi^ect of the bias resistors and Rh is essentially to 
provide an offset to the quasi- linear slope in Fig. 11. In the limit R^, ^ R^, 
the predicted values for 7_ would intersect (0, 0). The bias resistor, then, only 
provide a small correction by allowing a better fit to the measurements. 

Finally, there appears to be a maximum amount of damping where breathers 
can exist in the ladder. This maximum F occurs when 7+ coincides with 7_ 
which occurs at 0.2. This can be expected from the fact that the DB dis- 
cussed in this paper require hysteresis in the IV curves. 



5 Simulations at the experimental values 

5. 1 Introduction 

In section 2.1 we derived the standard model for the dynamics of the array. 
This model assumes the RCSJ for the junction dynamics, uniform bias condi- 
tion, and the effect of mesh self-inductances. We will use numerical simulations 
of the governing equations to provide for detailed analysis of the DB solutions. 
We numerically integrate Eq. (6) using a 4th order Runge-Kutta method. The 
initial condition is found via a similar procedure as in the experiments. We 
first bias the array near the depinning current. Then we increase the current 
of the middle vertical junction until it starts to rotate and finally we decrease 
this extra current to zero. Usually the resulting solution is periodic and we 
verify that it is stable by calculating Floquet multipliers. If the solution is not 
stable we integrate the equations of motion for a long time with an added 
small current noise source in every junction. In this way we perturb the solu- 
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Fig. 13. Simulation of type B breather with a subgap resistance modeled by Eq. 
(23). Prom bottom to top Qsg = 0.1,0.15,0.25,0.4, and 0.5. T = 0.18 and A = 0.05. 
When decreasing the current the type B breather destabilizes to a type A when 
Qsg = Q-'o- Inset shows the plot of Solid circles are from simulations and the line 
is from Eq. (16). 

tion and sample nearby trajectories in phase space. Usually the final result is 
a new periodic stable solution. 

Sometimes, and usually close to the destabilization points, when we integrate 
the equations we find aperiodic solutions that appear stable. We again add 
some noise to check the stability of the solution against fluctuations and use 
standard methods like Poincare section analysis or study of the Liapunov 
exponents of the system to gain more information about the behavior of the 
solution. 

5.2 Numerically integrated IV curves with a nonlinear resistance 

We will present numerical simulations of the dynamics of the ladder. In our 
samples h — 0.25 and the damping and the penetration depth vary with the 
temperature. Thus when the temperature varies from 4.2 to 9.2 K, F varies 
from 0.03 to 0.6 and A from 0.04 to 0.4. However, experimentally we only 
find breathers for T < 6.7K and that corresponds to 0.03 < F < 0.2 and 
0.04 < A < 0.05. We will then mostly present simulations at /i = 0.25, A = 
0.05, however the dynamics of the array is very rich and multiple transitions 
between different attractors occurs when changing the parameters. 

Fig. 13 shows simulated IV curves of single-site type B breather solutions. 
These simulations were done with a fixed F = 0.18 and A = 0.05. We have 
included a subgap resistance in this simulation in order to compare with the 
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experimental measurements more closely. We use the usual approach to extend 
the RCS J model where now R depends on the junction instantaneous voltage 
V{t). We define a conductance such that 



1 



ii\V(t)\<V, 



GiVit)) = 



RsaiT) 



(22) 



1 



otherwise 



Here Rsg{T) is taken to be only a function of temperature so for a given set of 
parameters it is constant. Our functional A/'((/9) now becomes (p+rg{v)(p+sm (p, 
where g{v) = R^/ Rsg{T) and V is calculated from R^. 

A simple approach is to model as a continuous hyperbolic tangent. We will 
use 



In our simulation, we take the value oi K = 100. Thereby Eq. (23) approaches 
a piecewise linear function with a conductance Qsg at v < 1 and 1 when v > 1. 

We find that the simulated curves are similar to the experimentally measured 
arrays. The inset compares /_ of the simulation to the prediction of the circuit 
model Eq. (16) when i?6 3> Rh- The deviation is due to inductance effects. 
Our simple circuit model neglects the effects of the inductances. We have also 
verified that /+ in the simulations is within 99% of the adjusted /+ of our 
model Eq. (18) and Eq. (21). We have found similar results for numerically 
computed type A breathers. 

There is little effect of the subgap resistance in our simulations besides chang- 
ing the shape of the IV curve. Instead of using a subgap resistance, we can 
redefine F so that it includes the effects of g. Then, V is calculated from Rgg 
and N'{(p) — (f + Tip + sin if as before. 

5.3 Breather instabilities and switching mechanism 

Another aspect of Fig. 13 is the resulting state of the array once the single- 
site type B breather becomes unstable. Some insight can be gained by first 
studying how the single-site breather destabilizes. Fig. 14(b) shows the dis- 
tribution of the Floquet multipliers at several values of the applied current 
close to /„ for the type B breather when g^g = 0.5 in Fig. 13. We use 
-^ext = 0.36, 0.35975, 0.3425. The voltage of the vertical junction V5 is 
approximately 0.55 in this small range of currents. Correspondingly, the 



9{V) = gsg + {I- ggsg) 



1 -tanhX(l -i;) 
2 



(23) 
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Fig. 14. Floquct multipliers of type A (a) and type B (b) periodic DB for decreasing 
currents above and at /_ of the simulation shown in Fig. 13 when ggg = 0.5. Figures 
(c) and (d) show as a function of the current the value of voltage (solid circles) 
and the modulus of the Floquet multipliers whenever the solution is periodic (open 
circles) . 



voltage of the top horizontal junction is V^t = V'5/2 — 0.275 V^. In our 
normalization, the fundamental frequency of the type B breather is then 
<^ = V'st/F = 1.53. Most multipliers lie on a circle of radius e"''^^''^/'^ Ki 0.83 
and this can be verified in Fig. 14(b). 

In Fig. 14(d) we decrease the current and show the value of the voltage and 
the modulus of the Floquet multipliers (only for periodic solutions). This is 
the typically observed bifurcation scenario for small A. It seems that for small 
A and underdamped junctions, this instability introduces more frequencies in 
the solution. When the periodic type B breather losses stability at / = 0.3425, 
the solution becomes a quasi-periodic type B breather similar to the one shown 
in Fig. 23. This quasi-periodic type B solution persists up to a current value 
of 0.33775 when the array jumps to a periodic type A solution. 

For large A, however, we usually observe a period- doubling bifurcation where 
a multiplier crosses the unit circle at —1, though the behavior also depends 
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on the damping. 

Our simulations in Fig. 13 sliow tliat tlie array sometimes switches to the su- 
perconducting state when the type B solution ceases to exist. The bifurcation 
scenario for the different subgap resistances is similar to Fig. 14(b). 

At Qsg = 0.5, on the other hand, when the single-site type B breather becomes 
unstable, the array switches to a single-site A breather. This type A breather 
exists for a range of currents. Fig. 14(a) shows the Floquet multipliers for a 
type A breather at current values: 0.27, 0.26975, . . . , 0.2495. The voltage of 
the vertical junction V5 is approximately 0.33 in this range of current. In 
our normalization, the fundamental frequency of the type A breather is then 
CO = V^/r = 1.81. Most multipliers he on a circle of radius e"^^^'^^'^ 0.85 
and this can be verified in Fig. 14(a). 

In Fig. 14(c) we decrease the current and show the value of the voltage and the 
modulus of the Floquet multipliers (only for periodic solutions). Below / = 
0.2495 the DB is unstable and the solution switches to the superconducting 
state. 

The selection rules between different m-site type B solutions shown by the 
experiments is an important feature that can not be explained in the frame- 
work of the DC circuit model and also is not well predicted by our simulations 
of the array. Once the breather solution reaches a critical value of the cur- 
rent the system chooses between many different attractors that have complex 
boundaries. In this process randomness and thermal fluctuations play a role 
as it is shown by the different switching patterns under similar experimental 
conditions. In our experiments we usually observe transitions from a m-site 
type B solution to a m + n-site type B. However sometimes we find jumps 
from type B to type A in the experiments (Fig. 9) The simulations (like in 
Fig. 13) show more frequently transitions from a type B solution to a type A 
solution and then, by decreasing the current, to the superconducting state. 

To try to understand this switching process of the cascade of m-site type B 
breathers, we have introduced more elaborated models to (i) include the bias 
circuit in our simulations, (ii) study the effect of external fields, (iii) introduce 
a full-inductance matrix formalism, (iv) take into account the nonlinear char- 
acter of the junction resistance, (v) include disorder randomizing the junctions 
critical currents and (vi) include thermal effects by adding a noise term to the 
junction currents. However, none of these simulations reliably predicts the ob- 
served cascade of m-site type B breathers found in our experiments as the 
current is decreased. 

The bias resistors can be added to the model by rewriting the KVL and 
KCL for the new circuit. We have also included an external magnetic field. 
However at small values of A the field is expelled from the array because of 
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strong inductive effects and does not affect tlie IV curves. Tfiis was also found 
experimentally in Fig. 10. 

It also relatively straightforward to use a full-inductance matrix. We find that 
the additions of extra couphng changes significantly the decay of the fields 
within the array. For instance, due to the non-local coupling the breathers are 
no longer exponentially localized, and the field decay is algebraic away from 
the breather. Nonetheless, the IV remains relatively unaltered. 

Finally, thermal effects may play an important role, especially if the attractors 
have complex boundaries. We use the standard Langevin approach and replace 
the resistor of the RCS J model by a noiseless resistor in parallel with a Johnson 
current noise source, 

CjVj + sin ^3 = ^xt + If (24) 

where < lf{t)I^ {t') >= {2kT/ Rj)6(t — t')6jk. This results in the usual current 
noise spectrum density Sj = 2kT/ Rj. We normalize our equations as in Sec. 2.1 
and N'{(pj) = (f j + Tipj + sin if j + and the spectrum of is Sj = 2kThjT /Ej 
where the Josephson energy Ej — {(^o/2tt)Icv and hj = ll/Icv 

Our dimensionless temperature is then T = kT/Ej. At 4.2 K our vertical 
junctions have lev = 345 /xA so that Ej = 8.2 x 10^ K. With this normalization 
4.2 K equals T ^ 5.1 x 10~^. Similarly, for the maximum temperature we 
observe breathers, 6.7 K, f 1.3 x 10~^. 

Our simulations were done aX h = 0.25 and A = 0.05. At these values of the 
parameters, we generally find in our simulations with Johnson noise that when 
a type B breather becomes unstable it jumps to a type A breather. Sometimes 
when a type B breather becomes unstable an m-site B breather will form as in 
our experiments, but this is much less common. However, we have found that 
at other values of the parameters (typically larger h and A) different switching 
scenarios occur. We have run several simulation that include thermal noise, the 
full inductance matrix, bias resistors, and possible stray fields. All these effects 
result in a very detailed model for the dynamics of the array and different 
breather solutions are found. Also, adding inhomogeneities, like a distribution 
of junction critical currents, or a magnetic field does not seem to change the 
fact that in the array simulations when most type B breather solutions become 
unstable a type A breather solution is formed. 

To investigate the experimentally observed cascade of m-site type B solu- 
tions we study a toy model that allows for only type-B breathers. The sim- 
plest model with this characteristic is a variation of the standard model 
where an up-down symmetry for the superconducting phases is assumed, thus 



29 



0.6 

mi 

c 

Fig. 15. Simulation of a cascade of 1-site, 4-site, and 6-site type B breathers in 
up-down symmetric ladder with T = 0.07 and A = 0.05. Each curve is the average 
of the indicated vertical junction. 

= — (^^(t) (this was also assumed in [26] and [40]). We remove the bottom 
branches of the full ladder, since their dynamic is by definition identical to the 
top horizontal junctions. Also, fluxoid quantization, Eq. (4), is modified since 

^] - - 2^5 + 27r/ + hf = 0. (25) 



This along with KCL yields a nonlinear coupled system that can be obtained 
from Eq. (6) defining f'j — and changing (^*- 2(^*- on the right hand sides. 

Fig. 15 shows a simulation of the up-down symmetric ladder with F = 0.07 and 
A = 0.05. The curves are the average voltages of each vertical junctions. We 
start at Igxt = 0.7 with a single rotating junction. As we decrease the current, 
the single-site breather becomes unstable at Jcxt = 0.32 and a 4-site breather 
is created. This solution persists until /gxt = 0.23 when a 6-site breather is 
formed and finally at I^xt = 0.15 the arrays switches back to the zero- volt age 
state. These current values are estimated well from our circuit model Eq. (16) 
and the discrepancies are due to the effect of the inductances. 

The up-down symmetric ladder only allows up-down symmetric solutions such 
as a type B breather. When the single-site breathers becomes unstable the ar- 
ray can no longer jump to an A type breather. It appears that this constraint is 
enough to allow the formation of m-site breather when the single-site breather 
becomes unstable. This toy model, thus allows for the study of the switching 
seen in the experiments, though it is not clear the physical relation of this 
model to the experiments. 
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6 Linear Analysis 



6.1 Resonances in a Ladder 



Before we embark on an analytical study of the different localized solutions 
in the ladder we first need to understand the basic linearized excitations that 
can occur. 

An important characteristic frequency of the ladder occurs when the frequency 
of a junction resonates with the lattice eigenmodes. To calculate the resonant 
frequencies, we linearize Eq. (6) around a solution. Every breather solution is 
approximately up-down symmetric far from the rotating junctions. Therefore 
we make the approximation = — (^^ and let (^*- = (pl+S(p*j and if'^ — ifg+Sip'j. 
The resulting linear equations are 



dip] + T6ip] + cos{ipl)6ip] = A(5(^J+i - 26cp] + 5ip)_^ + 25ip] - 25ip]_^) 



If the ladder is in the uniform whirhng state then the approximate solution is 
= and ^p^ = ujt + zj where z is the wavelength of the vortex train. To 
calculate the dispersion relation we let 6(f'j = e*e*(^-'+'^*) and 5(pj — e"e^^^^^'^\ 
We substitute in Eq. (26) and take the F = limit since our junctions are 
underdamped and it can be shown that for F < 1, the damping terms are only 
small corrections to the frequency. The cos((/9q)5(/9J results in terms of with 
coefficients of e*^"^* and and all other terms have coefficients of e*"^*. We only 
keep terms of e"^*. The resulting matrix equation is 



+ 1 + 



2A 
h 



2A(e-*^-l) -u;2^4Asin'(f) 






















(27) 



The solution to Eq. (27) is a;^ = F ± ^/F'^ - G where F = (1 + 2X/h + 
4Asin(2;/2)^)/2 and G — 4A sin(2;/2)^. Prom physical grounds we expect the 
wavelength to be well approximated hy z = 27r/, i.e. the average distribution 
of vortices in the array. The resulting resonant frequencies are important when 
studying properties of moving vortices in the ladder. Then, there is a traveling 
wave of vortices with density / that can resonate with the lattice modes. 

Our breathers solutions are clearly not uniform. So the above result, while 
instructive, is of limited value. We can think instead of a solution that is valid 
far from the localized breather core. Far from the breather core, the solution 
to first order is 



31 



(/7o = sin~^ iext 



(28) 



and this satisfies the ladder equations, Eq. (6). Since the core junctions in 
the breather are rotating, they will induce librations in all of the junctions 
in the array [26]. After substituting the first order solution into Eq. (26) and 
expanding the perturbations as plane waves, we are left with the following 
matrix that must have a zero determinant to allow for nontrivial solutions: 



uj"^ + iVuj + cos(7r/) + ^ 


A(e--l) 











2\{e-" - 1) _ 


l-iFo; + p + 4Asin^(|) 












(29) 



where p = yl — (icxt)^- We set / = as in the simulations and z represents 
the wavelength of the perturbations. This dispersion relation describes the 
linearized frequencies that can resonate with a breather. In the case where p 
is approximately 1 and F = the eigenvalues have a particular simple form 



/ 2A 

a;+ = W 1 + — + 4A sin^ {z/2) 
UJ_ = 1 



(30) 



we note that uj+ > U- and that c<j_ is associated with the LjC resonance of 
the ladder and uj^ with the LgC resonance. We will show that this is a good 
approximation for our simulations. 

A more generalized resonance condition can be calculated by using a harmonic 
balance technique. For the added complexity of the harmonic balance analy- 
sis, it only yields a result that is very similar to the linearized calculations. 
However, the harmonic balance approach can also be used to calculate the IV 
curve. 



6.2 Decay length 



To calculate a resonance frequency we substitute into (30) the appropriate 
wavelength. For instance, the resonance of largest frequency occurs when 
z — TT. However, whenever our frequency or wavelength falls outside of the 
dispersion relation then any linearized excitations must decay. In general we 
should let our wavelength be a complex number z — k + i^ where ^ represents 
the decay length. 



We first define sin^(2;/2) = x+iy and substitute into the determinant of matrix 
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Eq. (29). We expand terms and solve for x and y, 



4A(cj2 - 1) ATXlj 
4rXu; -4A(u;2 _ 





X 






y 





Here, the first row are the real components and the second row are the imag- 
inary ones. This equation can be solved for x and y, 



X — 



huj^ + [T^h - 2X - {2 + p)hy 



and 



4/iA[a;4 + (p^ - 2)uj^ + 1] 

[h(l + 2p) + 2A(1 + p) - V^hp - 2r2A]cj2 - [2A + h]p 
^ Ah\[u^ + (r2 - 2)u;2 + 1] 

Vuj{p-l) 



4A 2h[ou^ + {r-^ - 2)u^ + 1]' 

The solution simphfies drastically when p— I, 

hiu^ -2X-h 
x = - 



(32) 



(33) 



4A/i 



4A' 



(34) 



With these formula we can calculate the decay length. We use the fact that 
cos(A; + — cos k cosh ^ — i sin A; sinh ^. Then, 



X = K |sin^(2;/2)| = -(1 — cos A; cosh ^) 
y = $5 |sin^(2;/2)| = - sin /c sinh ^. 



(35) 



We can now solve for the k and ^. The solution is not simple in the general 
case but when y — 0, it simplifies considerably. When k is zero. 



^ = cosh ^ |1 — 2x\, 



(36) 



and the solutions decay exponentially. When ^ = we recover the normal 
modes of the system and the frequencies are given by the dispersion relation 
Eq. (30). 
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When y 7^ the solutions decay exponentiaUy for aU frequencies. This is due 
to the damping and is crucial for the existence of the resonant breather and 
chaotic breather solutions found in Sec. 7.2 and Sec. 7.3. 



Fortunately, our ladders are underdamped so the F — > limit is appropriate. 
Then, from Eq. (32) and Eq. (36), 



^ = cosh ^ 



-hu'^ + (2Xh + h{p + 1) + 2X)(jj^ - 2Xh - 2Xp - ph 
2Xh{cu^ - 1) 



(37) 



As p — > 1, this simplifies to 

^ = cosh"^ 



2X{h+l) + h 



2Xh 



2X 



(38) 



In the opposite limit of F — > oo, 1/^ approaches 0. 

Eq. (37) gives the result for the decay length of linearized excitations in the 
ladder when F = 0. Another limit where the equations simplify occurs when 



In the limit of a; —> we get 



cosh-' fl^llll^] , 

V 2Xh y ' 



(39) 



and it can be verified that this result is independent of F. 

We have taken the u; — > and the F ^ limits. There is one other important 
limit when inductances can be neglected (A oo). Then, with A ^ oo and 
p — > 1, the decay length is 

e = cosh-^ (40) 



and it can also be verified that this is independent of F. Also as A ^ 0, 1/^ 
approaches 0. 

The decay length given by Eq. (39) describes the decay of the DC flux in the 
array. For instance, using the parameters in Fig. 4 we calculate that 1/^ = 0.32 

and this agrees with the simulations. This results also implies that the decay 
is exponential and that this exponential localization has an upper bound. 
The decay length 1/.^ is always less than 1/ cosh~^{(/i + p)//i} if A 3> /i or 
l/cosh-^{(2A + p)/2A} if /i> A. 
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7 Numerical and analytical study of single-site DB solutions 

7. 1 Regions of existence of single-site DB solutions 

In section 5 we numerically simulated the behavior of the DB solutions found 
in the experiments. Such experiments were done at moderate to small values 
of the damping, anisotropy and penetration depth (F ~ 0.1, h ~ 0.25 and 
A ~ 0.05). It is also important to study the existence and behavior of the 
localized solutions at other values of the parameters and to estimate the critical 
values at which DB solutions destabilize. Varying the temperature we studied 
experimentally the behavior of the solutions at different values of the damping 
r. The results were presented in Fig. 11, where h ~ 0.25 and A ~ 0.05. For 
these values of the parameters good agreement was found with the theoretical 
predictions. 

The equations found in section 4 allow for a calculation of the IV curves and 
the maximum and minimum values of external currents supporting DB. For 
Ri, ^ Rh and single-site breather solutions we find 

N \ s ) R, 

for the IV curves and 

/+ _2h + s 
NIcv ~ 3h + s 

for the currents. 

Figs. 16 and 17 show the predictions given by the circuit model. The size of 
the existence regions decrease rapidly when the damping or the anisotropy 
increase. On the other hand, if the damping is small enough we should find 
localized solutions at large values of h. Also, the existence region is larger for 
type A solutions. 

This simple model, however, does not account for any dependence of the curves 
with the parameter A. This is an important limitation of the model and we 
have been unable to develop a more complete, yet still simplified, approach 
which incorporates A. We have confirmed in the numerical simulations that 
A affects our predictions in two important ways. First, it affects the value of 
the array retrapping current. The value used in our circuit models has been 
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Fig. 16. Prediction of Eqs. (42) for /-|- and /_ as a function of h for single-site type 
A (left) and type B (right) solutions and two values of the damping. Lightly hatched 
region corresponds to F = 0.08 and the densely one to F = 0.2 



1 

o 

I 

;o.5 



h=l / Type A 



h=0.25 



0. 



1 h=l / / TypeB 

-■-I / z 







0.25 

r 



s 0.5 iir 



01 



0.5 



h=0.25 



0.25 

r 



0.5 



Fig. 17. Prediction of Eqs. (42) for /-|- and /_ as a function of F for single-site type A 
(left) and type B (right) solutions and two values of the anisotropy. Lightly hatched 
region corresponds to h = 0.25 and the densely one to ^ = 1.0) 

calculated from a single junction and should be corrected by A in the case of the 
array. Indeed, some of the curves we will show below can be fitted assuming a 
simple linear dependence of the retrapping current with A. Second, as studied 
in the previous section, it governs the values of the voltage at which resonances 
between the breather and the normal modes of the array play an important 
role. Roughly speaking, the resonances split the diagrams in two different 
regions: The small and the large A regions. When A is small, the resonance 
frequency is smaller than the DB frequency, and when A is large the resonance 
frequency is larger. Thereby, complications of damped resonances between the 
DB and the lattice eigenmodes are avoided in these limits. 

Far from the resonance values the effect of A is a small correction to our IV 
curves. This is shown by the numerical simulations. See for instance Figs. 21 
and 26, where IV curves numerically integrated agree quite well with the 
predictions of the model, Eq. (41). 

We have done numerical simulations based on Eqs. (6) with / = in order 
to study the A dependence of the breather existence region. The results are 
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Fig. 18. Numerical calculation of the existence region of single-site DB when 
A = 0.04 and / = 0.6. Open circles correspond to type-A and solid circles to type-B 
solutions. Vertical lines correspond to cuts show in Fig. 20 and the asterisk to the 
experiments. 

presented in Figs. 18 through 20. In these diagrams we show the maximum 
and minimum values of the parameters for which a localized solution has 
been numerically found. As we will see, in some cases the characterization of 
the solutions inside the existence regions is quite complex in which several 
resonances and transitions between periodic and aperiodic states appear. In 
the figures we have also marked the values of the parameters at which the 
experiments were done, all far from these problematic areas and belonging to 
a region of the diagram where both type A and type B breathers are predicted 
to exist. 

The data were calculated by integrating the governing equations with a small 
quantity of noise. We start with a type B rotobreather and h ~ 0.001. As we 
increase h, type B solutions become unstable and the solution evolves to a type 
A rotobreather. As we further increase h this rotobreather becomes unstable 
and the system usually jumps to either a superconducting or a whirling state. 
To verify that our method is accurate, we have calculated Floquet multipliers 
for periodic rotobreather solutions and found results consistent with those 
shown. 

Fig. 18 shows the existence regions in the anisotropy versus damping plane 
when A = 0.04 and / = 0.6. When h is large the agreement with the predictions 
of Eq. (42) is good. There is an abrupt deviation of the curve for type B 
solutions at small values of h and a region at larger F where new type-B 
solutions appear. Vertical hues correspond F = 0.08 and F = 0.2, studied in 
Fig. 20. The asterisk corresponds to the experimental parameters. 

Fig. 19 shows the existence regions in the current versus A plane when F = 0.08 
and h — 0.25. These are the values of h and F in our experiments. We can see 
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Fig. 19. Numerical calculation of the existence region of single-site DB when 
r = 0.08 and h = 0.25. Open circles correspond to type-A and solid circles to 
type-B solutions. Horizontal lines correspond to the predictions of the circuit model 
and the asterisk to the experiments. 
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Fig. 20. Numerical calculation of the existence region of single-site DB when / = 0.6 
and r = 0.2 (left) and F = 0.08 (right). Open circles correspond to type-A and 
solid circles to type-B solutions. Horizontal lines correspond to the predictions of 
the circuit model and the asterisk to the experiments. 

that Eq. (42) gives a good estimation of the maximum and minimum values of 
the current for localized solutions, except for the case of the minimum current 
with a moderate A. In this case, resonances and other dynamical effects cause 
a substantial deviation from our simple model. 

Fig. 20 shows the diagram in the anisotropy versus A plane for / = 0.6 and 
r = 0.2 (left) and 0.08 (right). As expected from Fig. 16 localized solutions 
exist at larger values of h when F is smaller. Unexpectedly, at F = 0.2 and 
small values of A type-B solutions exist only for small values of h. This behavior 
is also described in Fig. 18. The asterisk in the F = 0.08 figure approximately 
correspond to the value of the parameter where our experiments were done. 



Clearly, the existence diagrams are complex and further research is necessary 
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Fig. 21. Simulated IV's for 9 x 1 ladder of type B breather as a function of A. 
h = 0.25, r = 0.1, / = and (a) A = 5, (b) A = 2, (c) A = 0.8, and (d) A = 0.5. 
The labels indicate different type B breather solutions. The vertical dashed lines 
are W-)-(A) and co- = 1 from Eq. (30). The horizontal dashed line is I_ from Eq. 
(16) when S> Rh and the diagonal dashed line is from Eq. (14) when Rj, » R^. 

to fully understand them. However, these diagrams show that DB solutions 
might be understood within our simple model at limiting values of small and 
large A. We will use these limits in the numerical and analytical analysis of 
type B and type A solutions in the following sections. 

7.2 Type B breathers 
7.2.1 Simulations 

Fig. 21 shows typical simulated IV's at different A's for single-site type B 
breathers. The horizontal dashed line is the minimum retrapping current ex- 
pected from Eq. (16) when i?5 ^ in the circuit model. The vertical dashed 
lines are the two resonant voltages of the ladder as approximated with Eq. 
(30) when z — (the largest lattice frequency). As can be seen in the IV's, 
a;+ gives a good approximation of the location of the resonant structure. Also, 
uj^ is always larger than c<j_. The diagonal dashed line is the expected voltage 
of the fifth vertical junction calculated from the DC circuit model, Eq. (14) 
with 3> Rh- The voltages of the fifth vertical junction and the fifth top 
horizontal junction are plotted and are usually related by V5 = 2V^t- Here, we 
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Fig. 22. Time evolution of the time derivative of the phase, v{t) = dip/dt, for the 
large A type B solution labeled Bla in Fig. 21(a) and the small A type B solution 
depicted as Bl in Fig. 21(d). We plot v{t) for the rotating vertical junction (5V) 
and two neighbor horizontal junctions (5T and 5B). Here, h = 0.25, F = 0.1 and 
/ = 0.8. 

have used the simple RCSJ model without a subgap resistance and we have 
set r = 0.1, the experimental value measured from the subgap resistance. 

Fig. 21(a) shows an IV when A = 5. It corresponds to the large A regime. 
At these parameter values the upper resonance is above any of the junction 
voltages. However, we still see some resonant behavior at / ~ 0.9. The solution 
when I — 0.8 is shown in Fig. 22(left). This type B breather is fully up-down 
symmetric in that </?* (t) = — </7^(i). Also, the averaged voltage of the horizontal 
junction is always half of the voltage of the vertical junction. Wc will label 
this type of solution as Bla. As wc decrease the current, the breather becomes 
unstable at / = 0.32 as predicted by the retrapping model, Eq. (16), and a 
type A breather forms which itself becomes unstable at / = 0.26. 

Fig. 21 (b) shows an IV when A = 2. We see that the V5 branch is now separated 
into two parts by the a;+ resonance while V^t is still below the resonance. The 
solution when all the voltages are below the resonance is still Bla. We will 
label as B 1/3 the breather solution where V^t is below and V5 is above the 
resonance. There is also a hysteresis loop that forms at the resonance that is 
not shown in the figure. We find that, as the Bla, the Bl/3 solution is also 
up-down symmetric, but in the case of the Bl j3 there exists a phase difference 
in between the velocity of the vertical and the horizontal junctions. 

When A = 0.8 we see the IV shown in Fig. 21(c). We find that we can interpret 

the f3 solution as a resonant type B breather. Below the resonance, there is a 
small remnant of the a breather. When the (3 solution becomes unstable but 
is still below the resonance, there is another type B breathers which we 
have labeled BI7. An unusual signature of this breather is that ^ ^V^t- 
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Fig. 23. Poincare sections of the third top junction (3T) for a BI7 solution. The 
phases are shown at times t = nr where r = 27r/V5 and h = 0.25, F = 0.1, 
A = 0.5, and / = 0.7. 

Further analysis shows that the 7 solution is an aperiodic type B breather. 
There is also a hysteresis loop at the resonance that, depending on the pa- 
rameters, might surround the aperiodic breather. This hysteresis would make 
it difficult to experimentally access this attractor by only using the applied 
current. 

We can continue to decrease A. Fig. 21(d) shows an IV when A = 0.5. The 
upper resonance has now divided the voltage of V5T into two branches. Below 
the resonance we get the solutions described above. We see that the /3 solution 
is getting "squeezed" by the a;+ resonance and the retrapping current. Indeed, 
there is a critical value of A where the type /3 solution disappears. 

Above the resonance, we find a solution labeled Bl. This small A regime is 
the same as in the experiments. The single-site breathers measured in our 
ladders are of this Bl type while the other breathers described above were 
only found numerically. In Fig. 22(right) we see that this Bl solution is not 
up-down symmetric although V^t = V5/2 = —V^b- 

We have done a Poincare section analysis of the aperiodic solution BI7 at the 
same parameter values as Fig. 21(d) In Fig. 23 we show the value of versus 
ip at times equal to to + nr where r = 27r/V5. The simplest periodic type B 
solutions have a period T = 2r since horizontal junction is half the voltage of 
the vertical junction. The solution shown in Fig. 23 seems to be quasiperiodic 
with two incommensurate frequencies, one of which seems to have a period 
equal to AT. 

There are many other type B breathers that were found numerically at other 
values of the parameters but are not discussed here. For instance, there is 
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Fig. 24. Poincare sections of the fourth and fifth vertical junctions for a type B 
chaotic solution. The phases are shown at times t = to + nr where r = 27r/V5 and 
h = 0.15, r = 0.2 and A = 0.2 and I = 0.6. (a) Shows the sections for the rotating 
vertical junction 5 and (b) for its first neighbor 4. 

a family of solutions which is not left-right symmetric. However, we shall 
focus our discussion to left-right symmetric solutions that have the above 
characteristics. 

At other values of the parameters, the IV curves show chaotic localized solu- 
tions. Fig. 24 shows Poincare section for vertical phases in the case of a chaotic 
localized type B solution. The values of the parameters for this solution are 
h = 0.15, X = 0.2, r = 0.2 and / = 0.6. Such chaotic region can be located in 
the central part of Fig. 20. 



7.2.2 Analysis 

In this section, we will use a harmonic balance technique to derive some ana- 
lytical descriptions of our DB. When studying periodic solutions it is almost 
always easier to work in Fourier space. Harmonic balance is a technique where 
the variables are decomposed into their Fourier components and substituted 
back into the governing equations. This creates a large set of coupled algebraic 
equations. If the governing equations are linear then each resulting algebraic 
equation is independent and the full system is easily solved. However, nonlin- 
earities tend to mix harmonic components. If the mixing effect is large then 
the resulting set of algebraic equations is usually intractable. However, in our 
underdamped Josephson arrays the capacitances act as filters that allow the 
transmission of only a few frequencies. Typically it is only one frequency. In 
this case, we can truncate the set of algebraic equations and a harmonic bal- 
ance technique can provide useful approximations. 

From simulations almost all breather solutions have left-right symmetry (and 
this is always true for average voltages) . We can then make a transformation 
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from the full ladder to what we call a half-ladder. A half-ladder is ladder where 
the breather is now on the first junction and we assume left-right symmetry. 

To make a mapping of the equations from a half-ladder to a full ladder we 
first consider placing a breather in junction 5 of a 9-j unction ladder. Current 
conservation at that node yields, 

r, =ii-il + P. (43) 

Mirror symmetry implies that — —I^ thereby current conservation becomes 

= 27* + r (44) 

If we have a half ladder, current conservation at the first node (labeled 5 for 
comparison) yields 

^ -Il + I' (45) 

To make a mapping between the circuit equations, we need to multiply /g 
by 2. Since P — hj\f{(p), this is equivalent to setting /ihaif ladder = 2/iiadder- 
However, wc want to maintain the same fiux pattern in the half-ladder so 
fiuxoid quantization must remain unchanged. This is simply done by setting 

Ahalf ladder 2Aiadder- 

Since the breather solution is highly localized we can, as a first approxima- 
tion, neglect every cell except the first one. The resulting reduced system of 
equations can be written as 



Af{cp') = -Af{cp') (46) 

hcAfiip')+Afiip')^ie.t (47) 

-hMi^')+Af{^n-ie.i (48) 

^^b^^i_^r^ ^U{v') (49) 



where he = 2/liadder and Ac = 2Aiadder- 

Most type B breathers have two voltages in the system corresponding to two 
frequencies. Our simulations indicate that the voltage of the horizontal rotat- 
ing junction can be approximately decomposed as 

ym = ±1 + a^'''^ COS (|t + + 6{*'^> cos (ujt + (50) 
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Similarly for the rotating vertical junction, 



/ 1 
V = UJ -\- a cos I — 




(51) 



For the librating junction 




(52) 



We integrate to get the phases and use exponential notation. Our phases are 
then 



where = (2a^/a;)e*^^ and = [If /Lo)e'^b with x = t,bj and r. By conven- 
tion we will take the real parts for the actual variables. Also, the integration 
constant has been taken to be k^*'^'^^ + 7r/2 for the rotating junctions. 

To substitute into our governing equations we must linearize the sine term. 
Substituting our phase ansatz into the sine term would yield an infinite Fourier 
series whose coefficients are Bessel functions of the amplitudes. In this sense 
the sine terms mix all of the harmonics. In the present case, our amplitudes 
are small and we can linearize cosx — 1 and sin a; = x. The sine for the top 
rotating phase is approximated as 




(53) 



sin(^* = K{e^(t*+'=*) 




'le^-k')} (54) 



and for the bottom phase 



sin if'' — K{e 




2 



+ ^e*(f (55) 



The sine of the vertical rotating phase is approximated as 



sin ^p■ = 3?{e' 




(56) 
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For the librating phase 

sin = 3?{-ie'^'' - i{A'-e'^' + cos T}. (57) 

We have neglected all terms of frequencies not equal to ±a; or ±a;/2. 

If we substitute our ansatz into our governing equations and expand the sine 
terms as indicated above, we transform the original 4 differential equations 
into a linear system of 20 equations with 21 unknowns. This extra degree of 
freedom is associated with the time translational invariance of the equations. 

In principle it is possible to solve the full algebraic system by fixing one of 
the unknowns, but here we will just estimate the amplitude oscillations of the 
breather. To make headway, we will use some trends found in the simulations 
to reduce the number of unknowns and take the limit F = 0. 

We first calculate and approximate the phases from the numerical simu- 
lations. From simulations, we find that as A ^ oo, the phase difference of the 
u>/2 harmonic between the waveforms of and sin ip^ is tt. We also find that 
the phase difference of the uj/2 harmonic between the waveforms of and 
is zero and between the waveforms of and (p^ is tt. These phase relations 
can be used to reduce the number of unknowns in the algebraic system. 

One solution to Eq. (46) obeys the up-down symmetry: — —B^ and — 
—A^. We can combine the uj/2 harmonic parts of Eq. (47), Eq. (48), and Eq. 
(49) and solve for \A\. This yields 

2h^\^-h^{uj/2f 
(u;/2)2(2Ae + 2/ieAe-/ie(a;/2)2) > 




Using the same technique to approximate phases using simulations for the uj 
harmonics, we find 



\B' 



1 



6^2(2 + 2//ie-CuV Ac) 



(59) 



with h\B*\ = \B^\. We can then find all of the amplitudes of the harmonics in 
terms of cu. 

As A decreases the breather enters a resonant regime. In this regime the phase 
relationships change from the above limit. We now find the phase difference of 
the ui/2 harmonic of the waveform of (/3* and sin is zero, the phase difference 
between the waveforms of 93* and if'' is tt, and between the waveforms of 99' and 
</7^ is also TT. Substituting the new ansatz in the governing equations results 
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1 2 3 4 5 6 7 




0.5 1 1.5 2 2.5 3 3.5 



A. 

Fig. 25. Solid circles arc simulated A* while dashed line is calculated A* from Eq. 
(58) for A large and Eq. (60) for A small. Solid lines is lo found in the simulation. 
Here T = 0.1 and h = 0.25. 

in the same equation as Eq. (58) but with an overall negative sign. The same 
occurs with the amphtude. In the resonant regime Eq. (59) just changes 
sign. 

In the A — >• 0, the equations simplify. From Eq. (49) we see that J\f{(p^) must 
tend to zero. This implies that 



This also puts strong constraints on A'^ and A . One possible solution is A'' = 
and A'' — 0. Then to satisfy Eq. (49), A* = A^ contrary to the expected up- 
down symmetry. However, this solution is consistent with Eq. (46) and it is 
the stable solution we find numerically for small A. 

Fig. 25 shows a summary of the analytical results and how they compare to the 
simulations. The top graph is a simulation for a cell while the bottom graphs 
shows results for our ladder. Both system have F = 0.1. For each system, we 
have excited a type B breather at iext — 0.5 and A = 10. We plot the value of 
1^4*1 and lu as we decrease A. 

As expected, both graphs are very similar when Aceii = 2Aiadder- The top solid 
line is from the simulated frequency of the rotating vertical junction. The 
bottom solid line is from the simulated frequency of the horizontal junction 
which is precisely half of the vertical frequency. The dashed line shows the 
a;+ resonance. Because there are two frequencies, the resonance divides the 
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c n 

Fig. 26. Simulated IV's for 9 x 1 ladder of type A breather as a function of A. Here, 
r = 0.1, / = and (a) A = 5, (b) A = 1, (c) A = 0.2, and (d) A = 0.02. The vertical 
dashed lines are w+(A) and = 1 from Eq. (30). The horizontal dashed line is /_ 
from Eq. (16) as 3> Rh and the diagonal dashed line is from Eq. (14) also as 

solution space into three regions that result in the three different solution 
found in Fig. 21 and studied analytically using harmonic balance. 

The solid circles in Fig. 25 are the simulated A* and the dashed lines are the 
analytical estimate. For small A we have the Bl solution whose amplitude is 
given by Eq. (60) and the predicted amphtude hes on top of the simulated 
values. For the resonance regime of solution Bl/3 the amplitude is given by 
Eq. (58), and for the large A regime the amplitude is also given by Eq. (58). 

For the simulations of a single cell, the harmonic balance gives a very good 
approximation. This is not too surprising, since the analysis was performed 
for a single cell. However, the bottom of Fig. 25 shows that the analysis is also 
valid for the full ladder as long as we transform A and h. 



7.3 Type A breathers 



Fig. 26 shows typical simulated IV's as a function of A for type A breathers. 
As in the previous discussion, the horizontal dashed line is the expected re- 
trapping current while the vertical dashed lines are the two resonant voltages 
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from Eq. (30). Here the voltage of the fifth vertical junction and the fifth 
horizontal junction are the same and this makes the analysis of the solutions 
much simpler. There are basically just two solutions. 

A harmonic balance technique can be used to calculate the amplitudes of the 
oscillating parts of the breather. We use a similar approach to Sec. 7.2.2 but 
approximate the phases with one frequency instead of two. Our phases are 
then 



- 3fJ + A;* + I - a*e^'"*| 
ip^ = ^ l^ut + + ^ - iA^e''^'^ 

(^'• = 3?{A;'^-aV'"*}. (61) 



When cu^ is above the junction voltages, as in Fig. 26(a) the solution is up- 
down symmetric with respect to the amplitudes: A* = —A'', and ~ A*. 
This means that all of the core junctions of the breather have a small oscil- 
lating amplitude. As A decreases the breather can resonates with the lattice 
eigenmodes as shown in Fig. 26(b). When A is so small that u;+ is below the 
breather frequencies we find a solution with A^ and ^ 0. 

The analysis of the amplitudes is straight forward and similar to what is done 
in Sec. 7.2.2. We substitute the ansatz in the governing equations and keep 
only DC and cu harmonics. We then use simulations to approximate some of 
the phase relations and thereby reduce the number of equations. 

In the A = limit, we have from Eq. (49) that A/'(v?*) = 0. If we set F = 0, 
this implies that 

1^*1 = (62) 



In the opposite limit, we can again use Eq. (49). We substitute the ansatz 
found in the simulations of A* = —A^ and approximate A^ as A*, then 



h/X 

3-hujyx 



(63) 



Fig. 27 shows a comparison the analytical results and how they compare to the 
simulations. The top graph is a simulation for a cell while the bottom graphs 
shows results for the ladder. Both system have F = 0.1. For each system, we 



48 




3 



ladder 



3,5 
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Fig. 27. Solid circles are simulated while dashed line is calculated from Eq. 
(63) for A large and Eq. (62) for A small. Solid line u) found in the simulation. Here 
r = 0.1 and h = 0.25. 

have excited a type A breather at iext = 0.4 and A = 10. We plot the value of 
I A* I and uj as we decrease A. 

As expected, both graphs are very similar when Aceii = 2Aiadder- The solid line 
is from the simulated frequency of the rotating vertical junction. The dashed 
line shows the uj+ resonance and we see that there are two types of solutions: 
Breathers with frequencies above a;+, and breathers with frequencies below 



The solid circles in Fig. 27 are the simulated A* and the dashed lines are the 
analytical estimate. For small A the amplitude is given by Eq. (62) and the 
predicted amplitude match the simulated values. For the large A regime the 
amplitude is given by Eq. (63) and the approximation breaks down at the 
resonance. In any case, we see that comparing a breather in one cell to the 
breather in a ladder works quite well. 



7.4 Vortex patterns in breathers 

The DB that we observe in the ladder can also be understood as bound vortex 
states. In Eq. (5) we define the vorticity rij in a cell as 



UJ+. 



ind 



(64) 
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where [</?] represents the phases modulus 2n, and fj^"^ = i^/27iX. As we will 
show, Uj changes in time differently for the three analyzed solutions. However, 
we want to stress that the magnetic field flux, which is due to the mesh 
currents in the ladder, librate around a mean value as shown in Fig. 4 and do 
not circulate. Only the fluxoid, or vorticity, which is the quantized quantity, 
moves from cell to cell. This behavior is different for the cases of type B or 
type A solutions and in the case of type B solutions for small or large values 
of A. 



7.4-1 Large A type B solutions 

In Fig. 28 we show a sequence of time snapshots of the ladder for a type B 
breather in the large A limit when / = 0. There is one vertical junction rotating 
and four horizontal rotating junctions. The voltage of the vertical junction is 
twice of the horizontal junction. For every full rotation of a horizontal junction, 
the vertical does two full rotations. The solid circle is a vortex and the "x" 
is an anti-vortex. In this large A case, the solution is completely up-down 
symmetric (see Fig. 22). Thus = — 

At the initial condition (a) there are no vortices in the ladder. As soon as the 
vertical junction goes over tt then a vortex- ant ivort ex pair is created as shown 
in (b) . This pair remains pinned in the ladder even though the applied current 
is large enough to depin isolated vortices. When the horizontal junctions go 
over TT then they create two vortices or antivortices as shown in (d) . The result 
is as if the vortex and antivortex pair change cells. This solutions remains until 
the vertical junction rotates through tt. Then, another vortex-antivortex pair 
is created that annihilates the pair of opposite polarity that was in the ladder. 
Now, there are no vortices in the ladder and the sequence repeats itself. The 
period double nature of the solution is evident. The vertical junction both 
creates a vortex-antivortex pair, and also annihilates the vortex-antivortex 
pair created by the horizontal junctions. 

7.4-2 Small A type B solutions 

Fig. 29 shows a type B breather in the small A limit at / = 0. This solution 
is similar to Fig. 28 but we no longer have up-down symmetry. Instead, the 
top and bottom rotating horizontal junctions have a n phase difference (see 
Fig. 22). We again start when the ladder has no vortices in (a). The vertical 
junction goes over tt and creates a vortex-antivortex pair. Then in (c) the 
bottom horizontal junctions go over vr and annihilate the vortex-antivortex 
pair. After some time, the vertical junction again creates a vortex-antivortex 
pair that disappears when the top horizontal junctions go over n and annihilate 
the vortex-antivortex pair. 
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(a) t=l, No vortices in system 
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(b) t=6. Vortex/anti-vortex pair created 
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(e) t=26 

XX *^ 



SX XX 

(f) t=30, No vortices in system 

Fig. 28. Vortex pattern of type Bla breather (defined in fig. 21) for large A at 

/ = o. 

In the small A limit vortex-vortex interactions are large. Therefore, vortices 
prefer not to enter through the top and bottom horizontal junctions at the 
same time. In order to minimize this repulsive interaction, vortices and anti- 
vortices enter the cells at different times when A is small and the up-down 
symmetry of the solution is broken. 
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(a) t=l, No vortices in system 



(b) t=6, Vortex/anti-vortex pair created 



Anti-vortex entered 



(c) t=13 



Vortex entered 



(d) t=22 



• \ X 



Kx 



(e) t=30, Vortex/anti-vortex pair created 
Anti-vortex entered Vortex entered 



/ 



T 

(f)t=34 

Fig. 29. Vortex pattern of type Bl breather for small A at / = 0. 
7.4-3 Type A solutions 

Fig. 30 shows a type A breather for / = 0. The solution in this case is eas- 
ier to interpret that the type B solutions since every junction rotates at the 
same voltage. First, the vertical junction goes over tt and creates a vortex- 
antivortex. The horizontal junctions then rotate over tt and annihilate the 



52 



(a) t=l, No vortices in system 



J L 



Vx 



(b) t=6, Vortex/anti-vortex pair created 
Anti- vortex entered Vortex entered 



If 



(c) t=13 

Fig. 30. Vortex pattern of type A breather at / = 0. 
vortex- antivortex pair as shown in (c). 

If we apply a field, then a phase difference develops between the horizontal 
junctions. For instance, in the appropriate field the right horizontal junction 
can go over tt first, thereby injecting a vortex to the right cell of the rotating 
vertical junction. Then, the vertical junction goes over tt creating a vortex- 
antivortex pair. The antivortex is created in the right cell and the vortex on 
the left. Since there was already a vortex on the right cell, the anti- vortex 
annihilates this vortex. The resulting state is as if the vortex moved from the 
right cell to the left cell. Then, the left horizontal junction goes over tt creating 
an anti-vortex in the left cell and thereby annihilating the vortex. This results 
in the ladder having no vortices. The full process can be interpreted as a 
rotating vortex around the common node of the rotating junctions. Therefore, 
type A breathers in a magnetic field are rotating vortices. 



8 Discussion and Conclusions 



Nonlinearity can cause localization in otherwise perfectly periodic systems. 

We have fabricated an experimental system to study these types of localized 
solutions called discrete breathers (DB). We have experimentally detected 
different one-site and multi-site DB states in a superconducting Josephson 
ladder network. DB can be excited and annihilated in a controlled way in the 
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ladder by using a local bias current. Then, by varying the external current 
and temperature we have studied the domain of existence and the instability 
mechanisms of these localized solutions. Experimentally we find two families of 
DB: type B and type A solutions. In the case of type B, four of the horizontal 
junctions in the array are in the resistive state, while for type A only two 
horizontal junctions are in the resistive state. We find that some of the type 
B and all the type A breathers do not obey the natural up-down symmetry of 
the ladder equations. The existence of type A solutions was predicted in [27] 
and recently has also been reported measurements showing a diversity (type 
B, type A and hybrid) of DB in an open J J ladder array [3]. 

We have developed a series of models and done theoretical analysis and sim- 
ulations to explain much of the behavior of the experimental data. At the 
simplest level, we have developed a DC circuit model for the system. In this 
model the junctions in the resistive state are studied as linear resistors and the 
junctions in the superconducting state as shorts. This model allows for evalu- 
ating the effect of the bias circuit and computing IV curves. We can estimate 
the minimum current for which both type A and B breathers can exist by us- 
ing this model and assuming that the minimum applied current is determined 
by a single junction retrapping mechanism. We can also predict the maximum 
current for which type B breathers can exist by using the junction's critical 
current. All these predictions agree with the experimental results. However, 
the type A maximum current found in our experiments is lower than predicted 
and is not yet understood. 

Using numerical simulations we have found the existence of many different 
DB attractors. We find periodic, quasiperiodic and chaotic type B solutions. 
Such solutions are identified with different regions in the numerically com- 
puted IV curves, and when we vary the parameters we find many different 
scenarios for the behavior of the solutions. The DB also excite waves in the 
ladder. We have calculated the frequencies of these small waves by using a 
linear analysis of the system equations. A resonance can occur when the DB 
frequency matches a frequency of the linear wave. We have numerically ver- 
ified that these resonances can cause instabilities in the localized solutions, 
but in most of the cases the resonances are not strong enough for destroying 
the localization of the breather. We have also performed a harmonic balance 
approximation which allows the characterization of the amplitudes of type A 
and type B breathers. 

Since DB do not exist for all the values of the model parameters, we have 
studied the regions of existence of localized solutions when these parameters 
vary. This study showed good agreement with the predictions of the simple 
model in many cases. In brief, DB exist at small values of the anisotropy, but 
the underdamped character of the junctions is crucial. In general, a smaller F 
results in a larger existence region. Thus, for small values of the damping DB 
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arc predicted to exist even at high values (> 1) of the anisotropy parameter. 
However, we find a complicated dependence of the DB instabilities with respect 
to A which is only partly explained by resonances between DB frequencies and 
the lattice linear waves. 

When varying A we find different physical behaviors for large and small values. 

At large values of A, type B solutions are up-down symmetric, and at small 
values of A this up-down symmetry is broken. In the intermediate A regions the 
resonances between the breather frequency and the linear modes are important 
and the behavior is much more complex. It is in this region where we found 
chaotic locahzed solutions. 

One important conclusion of these simulations at different values of A is the 
observation of resonances in the IV curves which do not destroy the DB solu- 
tion. This can be understood after the analysis presented in section 6.2. For the 
existence of DB, in the case of a forced and damped system the non-resonance 
condition is not necessary, because for any frequency the decay length is dif- 
ferent from zero. 

The induced magnetic fiuxes in the cells of the array librate around non- 
zero DC values and have different signs in opposite regions of the array with 
respect to the localization of the breather. To complete the physical picture of 
the dynamics of the DB we also analyzed the type B and type A solutions in 
terms of the vortex-antivortex behavior in the array. Briefiy, vortex-antivortex 
pairs are created in the center of the array and destroyed when new vortices 
and antivorticcs enter the array through the horizontal junctions. In the case 
of a type A solution with an applied external magnetic field, this vortex picture 
corresponds to the presence of a rotating vortex. 

Many of the aspects of the behavior of the solutions, including ranges of ex- 
istence, were understood using simple models. However transitions between 
different states when decreasing the external current is not fully understood. 
Although different behaviors have been found, a typical experimental result is 
that type B breathers become unstable through a cascade of switching events 
to multi-site type B solutions. To understand this phenomenon we have done 
more detailed simulations of the dynamics of the array. These simulations are 
based in the RCSJ model and we have numerically studied the effect of the 
bias circuit, the nonlinear resistance, full inductance matrix, thermal fluctu- 
ations and inhomogeneities in the junction parameters. Although sometimes 
a switching to a multi-site state was found, the typical numerical result for 
small A and h = 0.25 shows destabilization from type B to type A solutions. 
Another aspect is that all these transitions occurs close to the region of the 
gap voltage where it might be necessary to used models more sophisticated 
than the RCSJ one. 
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We have also studied a reduced model where up-down symmetry in the system 
is imposed by setting y^^j{t) = —^p'j{t). Then, transitions to multi-site states 
are typically found. However this model is unable to explain the existence of 
type A solutions and the transitions from type B to type A solutions that have 
been reported here and in [3] . This approach using a reduced model has been 
used in [26] and [40] to study the dynamics of the array. 

DB and vortices are two examples of localized excitations in a Josephson 
array. Future work will explore the interplay between DB and vortices (kinks 
or solitons) solutions which can coexist in the array [41]. Another interesting 
direction is a experimental study DB solutions in other geometries such as 
triangular arrays, and the excitation of breathers by AC driving currents. 

After this article was accepted we learnt of [42] which reports on numerical 
work of switching mechanisms in rotobreathers at /i = 0.44 and A ~ 0.21. 
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